Demultiplexing T cell activation experiment
10/7/2019
1 Introduction
Previously, we observed that time at room temperature (RT) prior to cryopreservation is a source of technical artifacts in blood-based single-cell RNA-seq (scRNA-seq) data. In this suproject, we aim to assess if culturing PBMCs after thawing them and before sequencing can remove the cold-shock signature.
2 Description of the data
To that end, we drew blood from four healthy donors (3 males and 1 female) and cryopreserved the samples after 0h (fresh), 8h or 24h at RT. Subsequently, we thawed the samples and either processed them immediately with 10X Chromium (day 0), or after 2 days of cell culture (day 2). To eliminate batch effects, detect doublets and reduce the library cost, we performed cell hashing. In this protocol, each condition (in our case time-points) is labeled with a specific hashtag oligonucleotide (HTO) that is crosslinked with an antibody. The antibodies bind to ubiquitous cell surface markers, and the HTO are sequenced alongside the single-cell gene expression libraries. We have the following libraries:
- Day 0: multiplexed with the two donors (male/female1 OR female2/female3), and the 3 time-points per donor (0, 8, 24h).
- Day 1/2: same as before but after 1 (BCLLATLAS_14 experiment) or 2 (BCLLATLAS_11) days in culture. Moreover, T cells were activated with anti-CD3 antibodies.
2.0.1 Objective
The objective of this notebook is to demultiplex the barcodes (cells) back to its original time-point. To achive that, we will follow the pipeline from Seurat.
3 Pre-processing
3.1 Load packages
library(scater)
library(SingleCellExperiment)
library(Seurat)
library(Matrix)
library(tidyverse)4 Demultiplex
# Load expression matrix, gene and cell metadata
libraries <- c(
"BCLLATLAS_11/Tcell_activation_day0",
"BCLLATLAS_11/Tcell_activation_day2",
"BCLLATLAS_14/Tcell_activation_day0_rep2",
"BCLLATLAS_14/Tcell_activation_day1_rep2"
)
t_act_list <- list()
for (lib in libraries) {
# Read the data
lib_path <- str_c("data/", lib, "/filtered_feature_bc_matrix/")
expression_matrix <- readMM(str_c(lib_path, "matrix.mtx.gz"))
barcodes <- read_csv(str_c(lib_path, "barcodes.tsv.gz"), col_names = FALSE)
colnames(barcodes) <- "barcode"
features <- read_tsv(str_c(lib_path, "features.tsv.gz"), col_names = FALSE)
colnames(features) <- c("ensembl", "symbol", "feature_type")
rownames(expression_matrix) <- features$symbol
colnames(expression_matrix) <- barcodes$barcode
# Separate HTO and RNA matrices
hto_ind <- which(str_detect(features$feature_type, "Antibody Capture"))
rna_ind <- which(str_detect(features$feature_type, "Gene Expression"))
t_act_hto <- expression_matrix[hto_ind, ]
t_act_rna <- expression_matrix[rna_ind, ]
# Setup Seurat object
t_act <- CreateSeuratObject(counts = t_act_rna)
# Normalize RNA data with log normalization
t_act <- NormalizeData(t_act)
# Find and scale variable features
t_act <- FindVariableFeatures(t_act, selection.method = "vst")
t_act <- ScaleData(t_act, features = VariableFeatures(t_act))
# Add HTO as an independent assay
t_act[["HTO"]] <- CreateAssayObject(counts = t_act_hto)
t_act <- NormalizeData(t_act, assay = "HTO", normalization.method = "CLR")
# Demultiplex
t_act <- HTODemux(t_act, assay = "HTO", positive.quantile = 0.99)
# Append to list of Seurat objects
t_act_list[[lib]] <- t_act
}
names(t_act_list) <- str_remove(names(t_act_list), "(BCLLATLAS_11/|BCLLATLAS_14/)")5 Visualize
We can visualize the results as ridge plots or heatmaps:
ridge_l <- map(t_act_list, function(t_act) {
Idents(t_act) <- "HTO_maxID"
RidgePlot(
t_act,
assay = "HTO",
features = rownames(t_act[["HTO"]])[1:6],
ncol = 2
)
})
heatmap_l <- map(t_act_list, function(t_act) {
HTOHeatmap(t_act, assay = "HTO", ncells = 5000)
})
# Save
# walk2(ridge_l, names(ridge_l), function(ridge, lib) {
# ggsave(
# filename = str_c("results/plots/", lib, "_hashtag_demux_ridge.pdf"),
# plot = ridge, height = 9,
# width = 16
# )
# })
# walk2(heatmap_l, names(heatmap_l), function(heat, lib) {
# pdf(file = str_c("results/plots/", lib, "_hashtag_demux_heatmap.pdf"), height = 4, width = 7)
# print(heat)
# dev.off()
# })ridge_l## $Tcell_activation_day0
##
## $Tcell_activation_day2
##
## $Tcell_activation_day0_rep2
##
## $Tcell_activation_day1_rep2
heatmap_l## $Tcell_activation_day0
##
## $Tcell_activation_day2
##
## $Tcell_activation_day0_rep2
##
## $Tcell_activation_day1_rep2
6 Save Seurat objects
# saveRDS(t_act_list, "results/R_objects/t_act_Seurat_list_demultiplexed.rds")7 Session Info
sessionInfo()## R version 3.6.1 (2019-07-05)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] forcats_0.5.0 stringr_1.4.0 dplyr_0.8.4 purrr_0.3.3 readr_1.3.1 tidyr_1.0.2 tibble_2.1.3 tidyverse_1.3.0 Matrix_1.2-18 Seurat_3.1.4 scater_1.14.6 ggplot2_3.3.0 SingleCellExperiment_1.8.0 SummarizedExperiment_1.16.1 DelayedArray_0.12.2 BiocParallel_1.20.1 matrixStats_0.55.0 Biobase_2.46.0 GenomicRanges_1.38.0 GenomeInfoDb_1.22.0 IRanges_2.20.2 S4Vectors_0.24.3 BiocGenerics_0.32.0 BiocStyle_2.14.4
##
## loaded via a namespace (and not attached):
## [1] readxl_1.3.1 backports_1.1.5 sn_1.5-5 plyr_1.8.6 igraph_1.2.4.2 lazyeval_0.2.2 splines_3.6.1 listenv_0.8.0 TH.data_1.0-10 digest_0.6.25 htmltools_0.4.0 magick_2.3 viridis_0.5.1 fansi_0.4.1 gdata_2.18.0 magrittr_1.5 cluster_2.1.0 ROCR_1.0-7 globals_0.12.5 modelr_0.1.6 RcppParallel_4.4.4 sandwich_2.5-1 colorspace_1.4-1 rvest_0.3.5 ggrepel_0.8.1 haven_2.2.0 xfun_0.12 crayon_1.3.4 RCurl_1.98-1.1 jsonlite_1.6.1 survival_3.1-8 zoo_1.8-7 ape_5.3 glue_1.3.1 gtable_0.3.0 zlibbioc_1.32.0 XVector_0.26.0 leiden_0.3.3 BiocSingular_1.2.2 future.apply_1.4.0 scales_1.1.0 mvtnorm_1.1-0 DBI_1.1.0 bibtex_0.4.2.2 Rcpp_1.0.3 metap_1.3 plotrix_3.7-7
## [48] viridisLite_0.3.0 reticulate_1.14 rsvd_1.0.3 tsne_0.1-3 htmlwidgets_1.5.1 httr_1.4.1 gplots_3.0.3 RColorBrewer_1.1-2 TFisher_0.2.0 ica_1.0-2 farver_2.0.3 pkgconfig_2.0.3 uwot_0.1.5 dbplyr_1.4.2 labeling_0.3 tidyselect_1.0.0 rlang_0.4.5 reshape2_1.4.3 cellranger_1.1.0 munsell_0.5.0 tools_3.6.1 cli_2.0.2 generics_0.0.2 broom_0.5.5 ggridges_0.5.2 evaluate_0.14 yaml_2.2.1 npsurv_0.4-0 fs_1.3.2 knitr_1.28 fitdistrplus_1.0-14 caTools_1.18.0 RANN_2.6.1 pbapply_1.4-2 future_1.16.0 nlme_3.1-145 xml2_1.2.2 rstudioapi_0.11 compiler_3.6.1 beeswarm_0.2.3 plotly_4.9.2 png_0.1-7 lsei_1.2-0 reprex_0.3.0 stringi_1.4.6 lattice_0.20-40 multtest_2.42.0
## [95] vctrs_0.2.3 mutoss_0.1-12 pillar_1.4.3 lifecycle_0.1.0 BiocManager_1.30.10 Rdpack_0.11-1 lmtest_0.9-37 RcppAnnoy_0.0.15 BiocNeighbors_1.4.2 data.table_1.12.8 cowplot_1.0.0 bitops_1.0-6 irlba_2.3.3 gbRd_0.4-11 patchwork_1.0.0 R6_2.4.1 bookdown_0.18 KernSmooth_2.23-16 gridExtra_2.3 vipor_0.4.5 codetools_0.2-16 MASS_7.3-51.5 gtools_3.8.1 assertthat_0.2.1 withr_2.1.2 sctransform_0.2.1 mnormt_1.5-6 multcomp_1.4-12 GenomeInfoDbData_1.2.2 hms_0.5.3 grid_3.6.1 rmarkdown_2.1 DelayedMatrixStats_1.8.0 Rtsne_0.15 lubridate_1.7.4 numDeriv_2016.8-1.1 ggbeeswarm_0.6.0